;; Figure the contours for 68.27%, 95.45%, and 99.73% confidence regions
FUNCTION DERIVE_CONTOUR_LEVELS, density, niter
  
  nlvs = long64(niter/10)+1
  level = dblarr(3)
  lvs = findgen(nlvs)/float(nlvs)*max(density)
  ret = dblarr(nlvs)              ; Set to zero
  totden = double(total(density)) ; For faster computing
  FOR pp=0LL,nlvs-1 DO BEGIN
     lind = where( density GE (lvs[pp])[0], nlind)
     area = (nlind EQ 0) ? 0.d : double(TOTAL(density[lind]))
     ret[pp] = area / totden
  ENDFOR
  lind1 = where(ret GE 0.682689492137d, nlevel1)
  lind2 = where(ret GE 0.954499736104d, nlevel2)
  lind3 = where(ret GE 0.997300203937d, nlevel3)
  
  level[2] = nlevel1 EQ 0 ? max(density) : max(lvs[lind1])
  level[1] = nlevel2 EQ 0 ? max(density) : max(lvs[lind2])
  level[0] = nlevel3 EQ 0 ? max(density) : max(lvs[lind3])
  
  RETURN,level
END
